model{
  ## observational model
  for (i in start[1]:end[1]) {
    y[i] ~ dnorm(mu1[i], nu1[i])
    mu1[i] <- alpha1[weekid[i]] + beta[house[i]]
    nu1[i] <- pow((tau[house[i]] * sqrt(mu1[i] * (1 - mu1[i]) / n.obs[i])), -2)
  }
  for (i in start[2]:end[2]) {
    y[i] ~ dnorm(mu2[i], nu2[i])
    mu2[i] <- alpha2[weekid[i]] + beta[house[i]]
    nu2[i] <- pow((tau[house[i]] * sqrt(mu2[i] * (1 - mu2[i]) / n.obs[i])), -2)
  }
  ## transition model
  alpha1[1] ~ dunif(0, 1)
  for (t in 2:n.weeks[1]) {
    alpha1[t] ~ dnorm(alpha1.star[t], omega.star[1])
    alpha1.star[t] <- psi[1] + phi[1] * alpha1[t - 1] + (t - 1) * xi[1] + 
      gamma[1] * x[t] + delta[1] * x[t] * alpha1[t - 1]
  }
  alpha2[1] ~ dunif(0, 1)
  for (t in 2:n.weeks[2]) {
    alpha2[t] ~ dnorm(alpha2.star[t], omega.star[2])
    alpha2.star[t] <- psi[2] + phi[2] * alpha2[t - 1] + (t - 1) * xi[2] + 
      gamma[2] * x[t + n.weeks[1]] + delta[2] * x[t + n.weeks[1]] * alpha2[t - 1]
  }
  ## prior distributions
  for (h in 1:10) {
    beta[h] <- beta.star[h]
  }
  beta[11] <- -beta.star[1] - beta.star[2] - beta.star[3] - beta.star[4] - 
    beta.star[5] - beta.star[6] - beta.star[7] - beta.star[8] - 
    beta.star[9] - beta.star[10]
  for (h in 1:10) {
    beta.star[h] ~ dnorm(0, 0.01)
  }
  for (p in 1:2) {
    psi[p] ~ dnorm(0, 0.01)
    phi[p] ~ dnorm(0, 0.01)
    xi[p] ~ dnorm(0, 0.01)
    gamma[p] ~ dnorm(0, 0.01)
    delta[p] ~ dnorm(0, 0.01)
    omega.star[p] <- pow(omega[p], -2)
    omega[p] ~ dunif(0, 100)
  }
  for (h in 1:11) {
    tau[h] ~ dunif(0, 100)
  }
}